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The concentration dependent diffusion problem is studied in two 
dimensions for a variety of cases relevant to the fabrication of narrow- 
channel mos transistors. A method for solving the partial differential 
equation is presented which allowed the solutions to be computed on 
an available minicomputer. This method enables analytic transfor- 
mations to be exploited which improve the accuracy of the solutions 
and increase the speed of computation. The simulations demonstrate 
that there are three principal nonlinear effects: (i) a large translation 
of the diffusion front, (ii) a marked steepening of the front itself, and 
(Hi) a very noticeable decrease in the ratio of the lateral to vertical 
diffusions. The ratio of the lateral to vertical diffusion as predicted 
by the model is compared with physical values experimentally deter- 
mined using a scanning electron microscope. 

I. INTRODUCTION 

In 1965, Kennedy and O'Brien analytically investigated the impurity 
atom distribution near the edge of a diffusion mask. 1 The geometric 
effects are significant, and the results of their study are routinely 
applied in semiconductor device design." The Kennedy and O'Brien 
model used the simple, linear, constant-coefficient, diffusion equation. 
In 1968, Hu and Schmidt introduced a concentration dependent dif- 
fusion model which included the effects of both the charged vacancy 
reaction and the impurity-induced electric fields. This model was 
investigated in one dimension, and it established that, in the regime in 
which semiconductors are actually fabricated, these nonlinear effects 
are quite significant for some impurities, such as arsenic. 



In the decade since these studies, there have been significant im- 
provements in fabrication technologies and lithographic techniques for 
both bipolar and mos transistors. Furthermore, the importance of mos 
transistors has also greatly increased. In this paper, we investigate the 
combination of nonlinear and geometric effects on the impurity atom 
distribution by studying a two-dimensional, concentration dependent 
diffusion model in a region containing the edge of the diffusion mask. 

In Section II, the nonlinear diffusion coefficient, including the effects 
of "autodoping," electric fields, and multiply ionized impurities, is 
derived in a form suitable for the calculation of the desired impurity 
profiles. The modeling parameters, the geometries, and the boundary 
conditions for the partial differential equations are also described. 
Section III outlines the numerical method used to solve the model. 
Additional details of the method are presented in the appendix. Several 
combinations of initial and boundary conditions are analyzed in this 
paper. Section IV analyzes the case where the impurity concentration 
in the window is held at a constant value throughout the diffusion, 
and Section V analyzes the case where a fixed quantity of impurity 
has been implanted and then diffused. This section also presents the 
simulation which is to be compared with experimental results. The 
simulations demonstrate that there are three principal nonlinear ef- 
fects: (i) a large translation of the diffusion front, (ii) a marked 
steepening of the front itself, and (Hi) a very noticeable decrease in 
the ratio of the lateral to vertical diffusions. 

To confirm the applicability of the physical model, a numerical 
solution of determined accuracy is compared with experimental mea- 
surements. Section VI contains a careful assessment of the accuracy of 
the numerical solution, while Section VII describes the experimental 
measurements of the lateral to vertical diffusion ratio. This experiment 
employed a scanning electron microscope (sem), and the resolution 
obtained is superior to that previously obtained by optical measure- 
ments. These measurements are compared with the predictions cal- 
culated in Section V. The comparison corroborates the validity of the 
concentration dependent diffusion model for arsenic and demonstrates 
its relevance for contemporary processing technology. Section VIII 
briefly discusses the application of the concentration dependent dif- 
fusion model to general process simulation. 

II. THE NONLINEAR, CONCENTRATION DEPENDENT DIFFUSION 
MODEL 

The diffusion can be modeled, as in Hu and Schmidt, 3 with a 
concentration dependent diffusion coefficient which incorporates the 
charged vacancy reaction and the self-induced electric field. For a 
noninteracting impurity, the transport of the impurity atoms in an 
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isotropic medium is described by a continuity equation for the impurity 
current, 

ir-- v - (J >- o) 

where N is the impurity concentration and J is the impurity current 
density. The current density is determined by diffusion and drift terms 
of the form 

J = -D N VN + ZhnNE, (2) 

where Dn is the diffusion coefficient, Z is the impurity charge, hn is the 
mobility, and E is the electric field. 

The self-induced electric field term can be approximated from the 
concentration, N, by assuming local space charge neutrality, 

n-p-ZN = 0. (3) 

Under these conditions, with the potential, ^, referred to the center of 
the band gap, solving 

n,exp (fl) " n ' exp (- f^J " ZN = 2n ' sinh (f£) - ZN = ° 



for \p yields 



kT 

yp = — S mh- ] (ZN/2n,), 
Q 



where n, is the intrinsic carrier concentration. Hence, 

E -v,.-*r ™ . (4 ) 

Q 2mJ(ZN/2m) 2 + 1 
Using Einstein's relation and this expression for E, (2) becomes 

J = -D N (l+ 1 Z ' N )VN. (5) 

V 2n,J(ZN/2n,y 2 + 1/ 

Following Ref. 3, we combine the electric field term in eq. (2) with the 
diffusion coefficient to obtain 

— = V.(D.vA(a)Va), (6) 

dt 



where 



h(a) = \ +-==== (7) 

>J(Za)' 2 + 1 
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and 

a = N/2n,. (8) 

Finally, we incorporate the following formulation of D N which was 
derived in Ref. 4 to account for the charged vacancy reaction, 

n -n {1 + (if) (9) 

Here Do is the phenomenological low concentration diffusion constant, 
f is the ratio of the electron or hole concentration to the intrinsic 
carrier concentration at the diffusion temperature ( /= n/n, for donors 
and f = p/rii for acceptors), and /? is a phenomenological coefficient 
discussed in detail in Ref. 4. By combining the equilibrium condition, 
np - n'f, with the assumption of local space charge neutrality, (3), we 
obtain the following approximation for /"in terms of the concentration 



f=\Z\a+ V(|Z|a) 2 +l. (10) 

Existing models of the type presented in this section have a long 
history of experimental validation. In his review paper, 5 Fair discusses 
in detail the cases of phosphorus, arsenic, and boron diffusion in 
silicon. In Ref. 3, Hu and Schmidt suggested that the correct values of 
fi were /? = 100 for donors and /? = 0.01 for acceptors. They also 
presented experimental evidence which tended to support these values. 
Subsequent experimental work, which incorporates more accurate 
measurements at high concentrations, has continued to corroborate 
($ - 100 for arsenic,' 57 but has shown that for boron ^ = 19 is the 
correct value. 689 For impurities such as phosphorus, which have mul- 
tiple vacancy reactions, this model is inappropriate. A model for 
phosphorus is discussed in detail in Ref. 10 and is beyond the scope of 
the present work. 

The assumptions upon which the present theoretical model is based 
are given in Ref. 4. Briefly, these assumptions are: (i) the density of 
the charged vacancies is much lower than the density of the impurity 
atoms, (ii) the charged vacancy population is in a state of quasi- 
equilibrium determined by the local equilibrium between the impurity 
atoms and the charged vacancies, and (Hi) the vacancy-impurity 
interactions are never more complex than pair interactions. These 
assumptions have been confirmed in the one-dimensional case by Fair 
and Tsai. 7 

In the present work, which is restricted to singly ionized impurities, 

Z=-\ and /? = 19 for acceptors, 

while 

Z = 1 and /? = 10" for donors. 
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Thus we want to solve 



— = V-(Z)(a)Va), 
dt 



(ID 



where the concentration dependent diffusion coefficient, D(a), is de- 
fined by eqs. (7) through (10), i.e., 



D(a) = Do 



1 + /3[a+ y/a 2 + lj 



1 + 



v/^TT 



:i2) 



Figures 1 and 2 present plots of D(a)/Do versus a for boron and 
arsenic respectively. Each figure displays the normalized diffusion 
coefficient for three ranges of a: 0.0 to 0.2, 0.0 to 2.0, and 0.0 to 20. 
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Fig. 1 — The nonlinear diffusion coefficient for acceptor impurities. 
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Fig. 2— The nonlinear diffusion coefficient for donor impurities. 

Although the behavior of the diffusion coefficient is mildly nonlinear 
for small values of a, it quickly approaches a straight line for a larger 
than 2. The limiting behavior of D (a) /Do for large a is 3.80a + 0.100 
for boron and 3.96a + 0.0198 for arsenic. This nearly linear behavior 
is consistent with the experimental results summarized in Ref. 9. It 
should be noted, however, that the approximations for both h{a) and 
f assume uncompensated material, whereas the experiments in Ref. 9 
are for heavily compensated material. 

The geometry of an idealized two-dimensional structure is presented 
in Fig. 3. It is a finite rectangle bounded by the lines x = xi, x = 0, 
y = yi and y = y/,. The mask, which runs along x = from v/ to 0, is 
assumed to be an impenetrable barrier for impurity atoms so that the 
flux across the interface is zero, i.e. da/dx = 0. The line y = y h is 
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Fig. 3 — The idealized structure. 



assumed to bisect the window so that the appropriate boundary 
condition is also that the flux across the edge is zero. Along the edges 
x = xi and y =yi the concentration will be held at a nominal background 
concentration, typically 5 to 7 orders of magnitude less than the initial 
peak concentration. The initial conditions and the boundary conditions 
in the window, x = and ^ y ^ y h , will depend on the specific 
problem being investigated. 

Two classes of diffusion problems are considered in this paper. The 
first class corresponds to the case where throughout the diffusion the 
concentration is held at a constant value, a () , in the window, x = and 
*£ y^yh . The corresponding initial condition is that throughout the 
interior of the semiconductor the concentration is set at a background 
level of 10 r, «o. This class of constant-surface-concentration diffusion 
problems is treated in Section IV. The second class corresponds to the 
case where an initial fixed quantity of impurity has been implanted, 
then diffused with a no-flux condition holding along the window. As 
an initial condition, we use an appropriately scaled Gaussian distri- 
bution corresponding to published results for actual implants. The 
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peak value of this distribution is denoted by a . This class of implant 
diffusion problems is treated in Section V. 

In each numerical solution which is described, the constants a , /?, 
and Do, as well as other relevant parameters, are specified and the 
time evolution is carried forward to a specified time, <„„„,. In discussing 
the solutions and their properties, it is frequently convenient to use 
the conventional dimensionless unit of diffusion length, where the 
dimensionless unit corresponds to physical length divided by 
v / 4Do<sto,.. This dimensionless unit is commonly used in diffusion 
problems which admit self-similar solutions. However, we should un- 
derscore the obvious fact that the problems we have just described do 
not admit self-similar solutions. This follows both from the bounded- 
ness of the domain and from the nonlinear nature of the differential 
equation. The implication, of course, is that each physical problem 
requires a separate numerical simulation, and in general this is true. 
However, in the case of constant-surface-concentration diffusion for 
wide mask devices it may be possible to develop a family of self-similar 
solutions parameterized by the surface concentration a . This is dis- 
cussed in Section VIII, which deals with the application of the nonlin- 
ear model to process simulation. 

III. THE NUMERICAL METHOD 

Equation (11) and its variations, which are derived later, are special 
instances of the following partial differential equation, 

f(t, X,y, U, U X , Uy, U,, Uxt, Uyl) = D X gl{t, X,y, U, U X , Uy, Ut, U xl , Uyl) 

+ D y g 2 (t, x, y, u, u x , Uy, u,, u xl , u yl ), (13a) 

with initial conditions 

u(0, x, y) = h(x, y) (13b) 

and boundary conditions 

a(t,x,y)n(x,y)-(g u g2) = b(t, x,y, u, u,), (13c) 

where the unknown function, u(t, x,y), is defined on the domain ktan 
*£ t ^ ktop, x, =s= x ss x h and yi < y «s y h , and where the vector, n{x, y), 
denotes the outward pointing normal along the boundary of the 
rectangle determined by xi, Xh,yi and yh. 

Given a well-posed problem of the form (13), the following numerical 
method is applicable. Represent the approximate solution, u(t, x, y), 
as a linear combination of tensor product B-splines with time-depen- 
dent coefficients 

m n 

u(t,x,y)= I I U Jll (t)B J (x)B k (y). 
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Using this representation, discretize (13) in space using Galerkin's 
method." This creates a large first-order system of nonlinear ordinary 
differential equations. The time evolution of this system can then be 
determined by any technique suitable for the solution of "stiff," implicit 
systems of differential equations. 1213 Our particular approach is to 
discretize in space with Galerkin's method applied to a second-order 
B-spline representation. An example of a second-order B-spline, usu- 
ally called a chapeau function, is shown in Fig. 4. In this figure the 
mesh points which define the chapeau function occur at integer coor- 
dinates (0, 0), (1, 0), (0, 1), etc., while the quadrature points used in 
Galerkin's method occur at (Vfe, Vfe), (— V6, Vfe), etc. The time evolution is 
then developed using extrapolation of the strongly A-stable implicit 
Euler rule. Additional details are presented in the appendix. 

By employing this general algorithm and by judiciously applying the 
chain rule for differentiation, we were able to develop our software so 
that the statement of the specific partial differential equation is 
separated from the general algorithm for solving it. There are two 
important consequences of this approach. First, it is extremely easy to 
investigate several models incorporating different physical effects. 
Second, it is easy to introduce analytic transformations of the original 
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Fig. 4 — A typical chapeau function. 
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problem which do not remove it from the canonical class defined by 
(13), but which may enhance the numerical accuracy and efficiency. 
Two such transformations are exploited in this paper. In Section 4.1 
we introduce a logarithmic transformation which substantially en- 
hances the accuracy of the solution. Then in Section 5.1 we introduce 
a Boltzmann-type, moving-coordinate system which enables us to 
obtain a solution in which the impurity dose is accurately conserved. 
All the simulations presented in this paper have been calculated on 
a coarse mesh of 11 X 21 points and on a fine mesh of 21 X 41 points. 
The meshes are defined as the cross product of 11 (21) points, uni- 
formly spaced between x, and 0, with 21 (41) points, uniformly spaced 
between y t and y h . The difference between the solutions on the coarse 
and fine meshes provides a reasonable estimate of the error for the 
coarse mesh. In several cases, asymptotic relations may be applied to 
this estimate to yield a good estimate for the error on the fine mesh. 
In addition to comparing the solutions on two different meshes, we 
have also applied checks against the linear problems, which have 
known analytic solutions, as well as against the nonlinear one-dimen- 
sional problems which can be solved numerically to very high accuracy. 
These checks on the numerical error are discussed in detail in Section 
VI. 

IV. CONSTANT-SURFACE-CONCENTRATION DIFFUSION 

In this section we discuss the numerical solution of the constant- 
surface- concentration diffusion problem, i.e., the case where the im- 
purity concentration in the window is held constant throughout the 
diffusion. Referring to Fig. 3, the initial conditions are that inside the 
semiconductor material (the rectangle determined by xi *£ x «= and 
yi *S y ^ yh), there is a background doping of 10" B ao, while along the 
window (x = 0, s£ y «s yh), the concentration is a () . The boundary 
conditions for this problem are the following. In the window the 
concentration is maintained at a,.. Along the sides x = x r and y = yi, 
the concentration is maintained at the background level, 10~ B ao. The 
mask is assumed to be an ideal barrier, so that the boundary condition 
along the mask (x = 0, y, =£ y *£ 0) is da/dx = 0. Finally, we assume the 
edge, y = y h , bisects the window so that the boundary condition along 
here is da/dy = 0. For each diffusion problem presented, we have set 
the constant Do to 1 and we have calculated the time evolution forward 
to time fcrtop = 1. 

4. 1 The logarithmic transformation 

The first point to observe is that the solution, which should always 
be positive, ranges over five orders of magnitude. This means that 
relative error— not absolute error— is the appropriate criterion. How- 
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ever, the Galerkin discretization in space tends to minimize absolute 
rather than relative error. This is also true of the vast majority of finite 
difference discretizations. One consequence is that an absolute error 
which is extremely small with respect to large values of the solution 
may be large enough to change the sign of small values. A relative 
error criterion is easily achieved if, instead of solving (11) directly, we 
derive and solve a partial differential equation for the logarithm of 
a(t, x,y) since, if a = exp(u) and a = exp(w + 5), then a a a(l + S) for 
small S. Let 

a(t,x,y) =e" u - x - y) 

and let 

D(u) = D(e"). 

Then (11) can be rewritten as 

u, - D(u)[u'i + uj] = V.(Z)(m)Vu). (14) 

Although this equation contains the additional nonlinear terms, 
u\ and uj, it can be solved more efficiently than (11), and the answers 
obtained are significantly more accurate in the sense of relative error. 

4.2 Linear, concentration independent diffusion 

The first case we consider is the case studied by Kennedy and 
O'Brien, namely, constant-surface-concentration diffusion for the lin- 
ear, concentration independent diffusion model on an idealized semi- 
infinite device. 11 Our interest in this case is twofold. First, this case 
offers a natural baseline against which to assess the significance of the 
nonlinear effects. Second, our numerical solution can be compared 
directly with the analytic solution provided by Kennedy and O'Brien. 
The accuracy of our discretization for this problem is then a rough 
estimate of the best accuracy achievable in the numerical solutions of 
the subsequent nonlinear problems. 

For this problem, the function D(ot) is simply the constant D , and 
eq. (11) simplifies to 

a, = V.(DoVa). (15) 

Let C(t, x, y) denote the solution developed by Kennedy and O'Brien 
in Ref. 1. This function satisfies equation (15) and the following 
constraints 

C(t, x,y)=0 for t = 0, x<0 

C(t,x,y) = l for t&Q, x = 0, 0*£y 

C x (t,x,y)=0 for t^O, x = 0,y<0. 

In order to use the Kennedy and O'Brien solution to check our 
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numerical technique, we pose a problem on the finite rectangle whose 
solution is C(t, x, y) + c . To do this, simply require that a x = along 
the mask and a(t, x,y) = C(t, x,y) + c„ along the remaining boundaries. 
Although this problem would form a mathematically precise test case, 
a somewhat more interesting test is to modify the boundary conditions 
to agree with those proposed for the nonlinear simulations, i.e., a(t, x, 
y) = Co along y = yi and x = x,, and a y = along y = yh. Inspection of 
C(t, x, y) shows that this problem is only a slight perturbation of the 
first problem if c = 10" 5 and the dimensions of the device, relative to 
tsio P , are greater than three diffusion lengths. 

The solution to the problem with Do = 1, t Mop = I, x, = y, = -5, and 
y h = 5 is presented in Fig. 5. This numerical solution is not only 
qualitatively correct, but the error in the concentration is less than 7.5 
percent, while the error in the log concentration is less than 2.8 percent. 
The precise definition of this error, which is relative to the maximum 
magnitude of the solution, is presented in Section VI. In that section 
it will also be pointed out that the larger errors are extremely localized 
and can be traced to two specific problems, the singularity at the edge 
of the mask and small displacements in the position of the diffusion 
front. 

4.3 Nonlinear, concentration dependent diffusion 

We now consider the nonlinear problem with the charged vacancy 
reaction and impurity-induced electric fields. Two cases of arsenic 
diffusion will be treated in this section. In the first case, a low- 
concentration case, the concentration in the window is held at a = 1, 
while in the second, high-concentration case, the concentration in the 
window is held at a ( > = 20. Because of the similarity of the diffusion 
coefficients, results for boron will be qualitatively the same. For both 
cases, Do and is...,, are set to 1. For the low-concentration case, the 
dimensions of the device are x/ = yi = -5 and y h = 5, while for the 
high-concentration case the dimensions are Xi = yi = —8.5 and yh = 8.5. 

Figures 6 and 7 present the results for a = 1 and a = 20, respectively. 
As might be expected, for the low-concentration case the diffusion is 
enhanced, but the differences from the linear case are modest. How- 
ever, for the high-concentration case the diffusion enhancement is very 
pronounced. The diffusion front has become extremely abrupt and has 
been translated from the region of 1 to 2 diffusion lengths from the 
mask edge in the linear case to a distance of 6 to 7 diffusion lengths in 
the high concentration case. 

The following rationale provides some qualitative insight into these 
nonlinear effects. When the concentration in the window is held 
constant, the diffusion coefficient near the window also remains con- 
stant and the solution in the immediate vicinity of the window behaves 
like the solution to a linear problem with an effective diffusion coeffi- 
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Fig. 5 — Constant surface concentration diffusion, linear model. 

cient of 3.96 Do a. However, this coefficient is falling linearly with the 
solution. In the critical region where a = 0.1 the solution becomes 
similar to that of a linear problem with a coefficient of Do. This causes 
the gradient of the solution to rapidly increase and generates an abrupt 
junction with a sharp corner. 

In summary, the combination of nonlinear and two-dimensional 
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Fig. 6— Constant surface concentration diffusion, nonlinear model with donor impu- 
rities and eta = 1. 

effects for the constant-surface-concentration diffusion implies that, as 
the surface concentration increases, the diffusion front is translated 
further into the device, the shape of the front becomes more abrupt, 
and the curvature in the corner region is increased. Moreover, these 
effects are modest for low impurity concentrations but extremely 
pronounced for high impurity concentrations. 
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V. DIFFUSION OF AN IMPLANTED IMPURITY 

In this section, we consider the case where an initial impurity dose 
has been implanted in the device just inside the window (x = 0, s£ y 
s£ y h ); see Fig. 3. The impurity is then diffused into the semiconductor. 
Along the mask {x = 0, y, s= y sS 0) and the edges x = x /t y = y,, and 
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Fig. 7 — Constant surface concentration diffusion, nonlinear model with donor impu- 
rities and on = 20. 
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y = y h , the boundary conditions are the same as for the constant- 
surface-concentration problem discussed in Section IV. However, the 
boundary condition in the window is now a no-flux condition, dot/dx 
= 0. In addition, we now have a nontrivial condition, which is the 
description of the impurity implant. 

In this section, we are concerned primarily with two simulations 
corresponding to a boron implant and an arsenic implant. The initial 
conditions for each implant will correspond to an appropriately scaled 
Gaussian whose mean and standard deviation are functions of the 
beam energy and the dose. We also assume that under the mask the 
implant decays quickly as a Gaussian in y. The case corresponding to 
a boron implant is somewhat hypothetical, but the case corresponding 
to an arsenic implant is carefully modeled after the key processing 
steps used to fabricate the junction studied in Section VII. 

5. 1 The Boltzmann coordinate transformation 

For this problem, it is useful to introduce yet another analytic 
transformation. A carefully selected moving coordinate system enables 
us to model the solution accurately throughout the transient evolution. 
This, in turn, improves our numerical conservation of the implanted 
dose. The natural transformation is 

y x t y * 

yJADot ylADot 

with the associated differential operators 

d 1 d d Id 



fo J4D t d £' & J4Dot d £' 
and 

dt~ dr 2t3£ 2t3£" 
This transforms eq. (14) into 



r* - I [to + to] - ^tt W + ?f] = v - 
2 4Do 



D(p) 

— — VI' 

4D 



(16) 



5.2 Linear, concentration independent diffusion 

Again, the first case we consider is the case studied by Kennedy and 
O'Brien, namely, instantaneous source diffusion for the linear, concen- 
tration independent diffusion model on an idealized semi-infinite de- 
vice. 1 Let C(t, x, y) denote the solution developed by Kennedy and 
O'Brien. This function satisfies (15) on the left half-plane with the 
boundary condition C x (t, x, y) = for t > 0, x = 0, and the initial 
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condition that it should yield an instantaneous unit dose per unit 
length for t = 0, y > 0. 

For our numerical problem, we consider eq. (15) with the initial 
condition 

aUstart, x,y) = CUsi„rt, x, y) + c 

and the boundary conditions 

«U» x, y) = Co for x = xi or y = yi 

a y (t,x,y) = Q for y = y h 

ot.x{t, x, y) = for x = 0, 

where kun = 10~ 2 and Co = 10~ 5 . Inspection of C(t, x, y) shows that if 
yh, relative to fetop, is larger than three diffusion lengths, than a{t, x, y) 
should be nearly identical to C(t, x, y) + Co. 

The solution of eq. (16) for Do = 1, t M()p = 1, x/ = yi = —5, and 
v/, = 5 is presented in Fig. 8. The key property of the Boltzman 
coordinate transformation is that the solution to (16) with the corre- 
sponding semi-infinite boundary conditions does not change its shape 
as time evolves. The amplitude of the solution simply decays as 
l/^JADat. Thus the solution to our finite problem should simply decay 
in the transformed coordinate system with no change of shape. This is 
indeed the case. The numerical accuracy proves to be limited only by 
the accuracy requested in the time evolution. For the example shown, 
the worst error is less than 0.4 percent. 

5.3 Nonlinear, concentration dependent diffusion 

We now consider two nonlinear problems, modeling diffusions of 
boron and arsenic implants. Referring to Fig. 3, the boundary condi- 
tions are that the concentration is held at the background doping along 
x = xi and y = yi. Along x = and y = yh, a no-flux condition must 
hold. If this condition is valid at t = t* {op , then for the transformed 
system we obtain 

"(t, £, = Co for f * & or £ = £/ 

v s (r, f,£)=0 for f = 

Mt,£|)=0 for £ = £,,, 

where t 2= Tsi lir t . 

For the boron simulation we consider a hypothetical device with 
Xi = y, = —2.0 |im and yh = 2.0 /xm. The diffusion temperature is taken 
to be 1400 °K. The implant is taken to be the familiar Gaussian with a 
spread of 0.04 /mi and a peak concentration of ao = 18.45, where the 
concentration is normalized to the intrinsic carrier concentration. The 
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Fig. 8 — Implant diffusion, linear model. 

background doping is set seven orders of magnitude below this peak 
concentration. The length of the diffusion is 900 seconds, and the value 
of kiart is set at 65.0 seconds. The choice of fci Hr , is dictated by the 
resulting value of the low-concentration diffusion constant, which in 
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this case is D = 3.37 X 10~ 13 cnr/s, and the requirement that the 
boundary conditions will be valid on the window defined by f/ £/, and 
£h ■ The result of the boron simulation is shown in Fig. 9. 

For the arsenic simulation, the dimensions of the device are taken 
to be xi = y t = -0.9 /xm and y h = 0.9 fim. The diffusion temperature is 










Fig. 9 — Implant diffusion, nonlinear model for boron. 
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1373°K. The implant is assumed to have a spread of 0.017 p.m and a 
peak concentration of 114. Again, the background doping is set seven 
orders of magnitude below the peak concentration and the length of 
the diffusion is 900 seconds. In this case, Do = 2.03 X 10" N cnr/s, and 
f slar1 is set to 67.2 seconds. The result of the arsenic simulation is shown 

in Fig. 10. 

As the peak concentration drops, the nonlinear model becomes 
increasingly linear. However, the solutions shown in Figs. 9 and 10 
clearly demonstrate that for contemporary fabrication the combination 
of peak concentration and diffusion time cause the solutions to exhibit 
several nonlinear effects. In particular, the shape of the diffusion front, 
the curvature of the corner, and the ratio of lateral-to- vertical diffusion 
distances are at variance with the linear predictions even for boron. 

VI. NUMERICAL ACCURACY 

A numerical simulation should be subjected to the same careful 
scrutiny as a physical experiment to establish error bounds and insure 
that the results are free from any systematic bias. This is particularly 
important in the current study where we are employing a new algo- 
rithm on a relatively small computer to investigate the two-dimen- 
sional effects predicted by a nonlinear model. In this section, we 
discuss the checks we have applied to this study and present various 
estimates for the accuracy of the numerical results. 

First, we establish some meaningful but concise error descriptions. 
Let a(t, x, y) denote the concentration and u(t, x, y) the log concen- 
tration, exp(u(t, x, v)) = a(t, x, y). Let ||a|| and ||u|| denote the 
maximum magnitude of each function over the rectangle x, *£ x =s= 
and yi^y^yi,. If a and u are computed solutions, define e(t, x, y) = 
| 6t{t, x, y)-a(t, x, y) | and 8(t, x,y) = \ u(t, x, y) - u{t, x, v) |. Then two 
concise error measurements are || e || / j| a || and || 8 1| /|| u || , where || e || and 
|| 6 1| denote the maximum magnitude of these functions over the 
rectangle. It is important to realize that these conventional criteria, 
although mathematically precise, represent the worst-case error and 
may be unduly pessimistic or even misleading. 

We have applied three types of checks to the simulations. First, we 
have computed the numerical solution to linear problems for which we 
already know an excellent analytic approximation. This serves first as 
a consistency check and second as a gauge of the best achievable 
accuracy for the nonlinear results. Second, we have computed each 
solution on a coarse mesh of 11 x 21 points and on a fine mesh of 21 
X 41 points. Four-thirds of the difference between these two solutions 
yields a good estimate on the error in the coarse mesh solution. With 
some exceptions which are discussed later, the error in the fine mesh 
solution should be one-third the difference. These estimates of the 
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Fig. 10 — Implant diffusion, nonlinear model for arsenic. 

error follow from the fact that we are applying a second-order discre- 
tization in space. Finally, the solutions along the line y = y h should 
correspond to the solution of a one-dimensional problem, provided, of 
course, that the dimensions of the device are sufficiently large. Our 
third check is to compute the solutions to these one-dimensional 
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problems to very high accuracy. We use a very well-tested one-dimen- 
sional pde package, post, 14 and meshes of 11, 21, 41, and 81 points. 
These highly accurate solutions are then compared with the two- 
dimensional solutions to provide a check along this boundary. This 
last approach yields a consistency check both on the numerical soft- 
ware and on the dimensions of the device in that the one-dimensional 
numerical solution for 21 mesh points should be identical to the 
corresponding two-dimensional solution. In every case, the one- and 
two-dimensional solutions for the same mesh do, in fact, agree to the 
precision requested in the time evolution. 

6. 1 Errors in the implant diffusions 

Since the error analysis is simpler in the case of the implant 
simulations, we address them first. Figure 11 graphs the four one- 
dimensional solutions for boron. A number of interesting facts can be 
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Fig. 11— One-dimensional solutions, nonlinear model for boron. 
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Fig. 12— Relative difference in the concentration, nonlinear model for boron. 



observed. First, the solutions are rapidly converging as the mesh is 
refined. In fact, closer inspection shows that the asymptotic h 2 rate of 
convergence is being achieved, that is, the error is being decreased by 
a factor of one-fourth each time the mesh spacing is halved. Second, 
the 21 -point solution, while still discernibly in error, is quite good. The 
maximum error in the concentration is only 9.6 percent, while the 
maximum error in the log concentration is only 8.4 percent. Third, we 
see that for boron the chief contribution to the error is an underesti- 
mate of the amplitude. This is probably due to the fact that the coarser 
meshes do not conserve the dose well. We also note that the maximum 
error in the log concentration occurs in the steepest part of the front 
where a small change in the x coordinate introduces a very large 



CONCENTRATION DEPENDENT DIFFUSION 23 



change in the value of the concentration. In fact, it is probably far 
more relevant to remark that the maximum error in the position of 
the front is less than 0.06 /*m, a relative error of about 4.0 percent if 
the front is assumed to be located 1.5 /xm inside the window. 

Figures 12 and 13 are plots of the relative differences between the 
coarse (11 x 21) and fine (21 X 41) mesh, two-dimensional solutions. 
Both figures show that the maximum error is quite localized and that 
it should be well estimated by the one-dimensional errors. Conse- 
quently, we can safely estimate the maximum errors in the fine mesh 
solution by applying the h 2 relation, i.e. $&». ~ (Wcoaree — wnne)/3. The 
error estimates for boron are summarized in Table I. 

Because of the increased nonlinearity of the donor model as well as 
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Fig. 13— Relative difference in the log concentration, nonlinear model for boron. 
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Table I — Error es timates for the implant simulations 

Error Estimate Boron Implant Arsenic I mplant 

One-Dimensional 

« 11/ II « II 
Two-Dimensional 

\U\\/\\o\\ 
Frontal Displacement 
Absolute error 
Relative error 



8.4% 
9.6% 


2.4% 
2.2% 


5.6% 
5.5% 


32.6% 
33.0% 


0.06 fim 
4.0% 


0.018 (im 

2.7%. 
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Fig. 14 — One-dimensional solutions, nonlinear model for arsenic. 

the much higher initial peak concentration, the case of the arsenic 
implant is considerably more difficult numerically than boron. Figure 
14 shows the four one-dimensional solutions for this case. Two facts 
leap out. First, the 11-point solution is extremely inaccurate, while the 
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21 -point solution is in fact more accurate than the corresponding 
solution for boron. The maximum estimated errors in the concentra- 
tion and log concentration are 2.2 percent and 2.4 percent, respectively. 
Second, the front is somewhat steeper which, in turn, amplifies the 
sensitivity of the error to small changes in the position of the front. 
The front is displaced by about 0.018 /xm, a relative error of about 2.7 
percent if the front is assumed to be located 0.67 /xm inside the window. 
In spite of the relatively good accuracy, the errors indicate that the 
asymptotic h 2 rate of convergence has not yet been achieved. 

Figures 15 and 16 are plots of the relative differences between the 
coarse (11 X 21) and fine (21 x 41) mesh, two-dimensional solutions. 
Both plots indicate that the dominant error is caused by the displace- 
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Fig. 15— Relative difference in the concentration, nonlinear model for boron. 
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Fig. 16— Relative difference in the log concentration, nonlinear model for arsenic. 

ment of the front. Figure 16 also indicates that the error in the log 
concentration should be well estimated by the one-dimensional anal- 
ysis. 

Table I summarizes the error analysis for the implant study. The 
one-dimensional and frontal displacement error estimates are taken 
directly from Figs. 1 1 and 14. The two-dimensional error estimates are 
simply the maximum differences shown in Figs. 12, 13, 15, and 16, 
multiplied by the asymptotic correction factor, %. In the case of boron, 
the one-dimensional analysis indicates that, although the errors are 
beginning to display the proper asymptotic behavior, the two-dimen- 
sional estimate is overoptimistic. In the case of arsenic on the other 
hand, the one-dimensional analysis shows that, although the errors are 
not yet behaving asymptotically, the two-dimensional estimates are 
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grossly pessimistic. Clearly, the poor placement of the front in the 11 
x 21 point mesh completely dominates the error. We conclude that 
the 21 X 41 point mesh is quite adequate for engineering work. We 
also note that merely reporting the maximum error is too simplistic. 

6.2 Errors in the constant-surface-concentration diffusions 

In the case of constant-surface-concentration diffusion, the concen- 
tration in the window is pinned throughout the diffusion. This means 
that the dominant error in the one-dimensional region is simply the 
frontal displacement. However, this also means that the true solution 
is singular at the edge of the mask {x = 0, y = 0). Such boundary 
singularities invalidate the h 2 asymptotic property of the numerical 
method we are employing (as well as virtually any other numerical 
method that does not deal with the singularity directly). Fortunately, 
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Fig. 17— One-dimensional solutions, constant surface concentration diffusion, a« - 1. 
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Fig. 18 — One-dimensional solutions, constant surface concentration diffusion, on = 20. 



the effect of such a singularity is generally restricted to the immediate 
vicinity of the singularity. 

Figures 17 and 18 present the four one-dimensional solutions for the 
low and high concentration cases, respectively. The frontal displace- 
ment problem is particularly acute in the high concentration case. The 
error in the position of the front is about 11.6 percent. While this may 
well be acceptable as far as locating the front, the fact that the front 
is so steep means that the value of 8/\\ u \\ in this region is essentially 
meaningless. 

Another way of describing the error in a problem of this nature is to 
apply a type of backward error analysis. It may be more informative 
to compare the solution with a more accurate solution evolved forward 
to a different time. For each solution of the four constant-surface- 
concentration cases, we have picked a mesh point in the diffusion front 
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and have then required that the 11-, 41-, and 81-point solutions take 
on the same value at this position as the 21-point solution. Figure 19 
presents the results for this one-dimensional study in the high concen- 
tration case. From this point of view, the error in the stopping time is 
18.1 percent, while the maximum error in the concentration is only 1.0 
percent, and the maximum error in the log concentration is 2.3 percent. 
We now consider the effect of the mask edge singularity. That the 
mask singularity has an effect is clearly evidenced from the small 
irregularities near the mask edge in Figures 5 to 7. Figure 20 presents 
the relative error in the concentration between the numerical solution 
to the linear problem and the analytic solution. Two things are 
apparent. First, the major contribution to the error is the mask 
singularity, and second, its effects are quite localized. This second 
observation is important if we are to place any confidence in the 
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Fig. 19— One-dimensional solutions, L varying, constant surface concentration 

diffusion a„ = 20. 
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Fig. 20— Relative error in the concentration, constant surface concentration diffusion 
linear model. 



prediction of the position, y ( (y, < y, ^ 0), of the "lateral out-diffusion 
front." Figure 21 presents the relative error in the log concentration 
for the linear problem. Figure 21 indicates that, for the nonlinear 
model, the dominant error in the log concentration will be due to the 
frontal displacement. Comparisons of the course and fine mesh solu- 
tions yield plots that are qualitatively the same as the last two figures. 
The error in the concentration is dominated by the mask edge singu- 
larity, and the error in the log concentration is dominated by the 
frontal displacement. 

The error estimates for the constant-surface-concentration diffu- 
sions are summarized in Table II. The absolute error in the location of 
the diffusion front is given in diffusion lengths. The relative time error 
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Fig. 21— Relative error in the log concentration, constant surface concentration 
diffusion, linear model. 



Table II — Error estimates for the constant-surface- 
concentration simulations 



Error Estimate 


an = 1 


a„ = 20 


One-Dimensional 






1 


S 11/ II Mil 
f / a\\ 


1.2% 


2.3% 




0.3% 


1.0% 


ht/t Ht „ v 


5.9% 


22.2% 


Two-Dimensional 






If KIWI 


3.6% 


34.1% 


12.8% 


9.6% 


M/t slup 


9.7% 


37.2% 


Frontal Displacement 






Absolute error 


0.08 


0.85 


Relative error 


3.6% 


11.6% 
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estimates for the two-dimensional cases are determined by evolving 
the coarse mesh solution forward until it interpolates the fine mesh 
solution at a prescribed point in the diffusion front. The one-dimen- 
sional error estimates are simply the maximum differences between 
the 21-point and 81 -point solutions, and the two-dimensional estimates 
are simply the maximum differences between the 11 X 21 point and 
the 21 X 41 point solutions. Obviously, the two-dimensional estimates 
are unduly pessimistic. However, no asymptotic relations appear to be 
applicable because of the difference in stopping times as well as the 
presence of the mask edge singularity. Again it seems reasonable to 
conclude that the accuracy obtained by the 21 X 41 point mesh is 
adequate for engineering work. 

The analysis presented in this section would appear to conclusively 
establish four facts: (i) the algorithm outlined in Section III is effective, 
(ii) the various diffusion simulations are free from systematic errors, 
(Hi) for acceptor and donor impurities in low concentrations, the 
simulators can obtain answers accurate to a few percent, and (iu) for 
high concentrations of impurities, the simulators can obtain answers 
which are qualitatively correct and probably suitable for engineering 
design. At the same time, the fact that for high concentrations the 
errors were on the order of 10 to 30 percent clearly indicates that more 
accurate results in this regime will require more sophisticated tech- 
niques or larger computers. 

We have indicated two analytic transformations of the partial dif- 
ferential equation which have made it more tractable numerically. 
Several other ideas have also suggested themselves and should be 
explored some time in the future. First, boundary singularities, such as 
the one caused by the mask edge, are well-known problems, and for 
simple linear problems there are two standard remedies, either explic- 
itly model the analytic singularity or locally refine the mesh in the 
vicinity of the singularity. For nonlinear problems, the same ideas 
should apply and we have briefly explored both possibilities, but we 
have not as yet achieved any notable success. Second, the underlying 
mesh used in the calculations presented in this paper consists of points 
uniformly distributed along the x and y axes. Clearly, a graded mesh 
could further enhance the accuracy of the solution. However, while a 
graded mesh will not adversely affect the asymptotic behavior of the 
errors, care must be exercised in grading the mesh to insure that the 
conditioning of the Galerkin matrix is not unduly degraded. Third, in 
view of the generality of the underlying differential equation, a natural 
idea would be to explore curvilinear coordinate systems. The level 
curves of the solutions immediately suggest a system based on the 
conformal map, w = (z + 2 -1 )'/-, which gives rise to an elliptic- 
hyperbolic coordinate system. 
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VII. PROFILE MEASUREMENTS 

The calculations of concentration dependent diffusion profiles dis- 
cussed here depend, for their validity, on the applicability of the 
underlying physical theory developed by Hu. 7 Experimental measure- 
ments of the vertical diffusion have verified the one-dimensional 
version of the model:' Verification of the two-dimensional model 
requires that the lateral diffusion effects also be examined. This section 
presents the results of such an experiment. 

Measurement of lateral diffusion effects by angle lapping is impos- 
sible since the vertical extent of the doping profiles precludes magni- 
fication by angle lapping. The resolution of optical microscopes is 
insufficient to provide accurate measurements of the small lateral 
distances involved, ~0.5 /un. The resolution problem is solved by 
making measurements of required distances using a scanning electron 
microscope (sem). The remaining difficulty arises from exposing the 
junctions to be measured and preparing samples in a way which will 
allow the junctions to be detected. 

Two methods of junction exposure were attempted. Conventional 
angle lapping of the junction has been previously used to expose 
junctions. 15 This technique results in smooth, easily cleaned surfaces, 
but control of the specific location of the surface cut is difficult. A 
second method of laser scribing to ±25 /un on either side of the junction 
and cleaving provides accurate ±25-jum surface cut location, but leaves 
a surface which can contain numerous surface cracks which trap 
charge and obscure the surface during sem studies. Both methods 
allow operation of the device with substantially the same electrical 
characteristics as before the sample preparation. 

Two kinds of imaging techniques were used to make the junction 
position visible. MacDonald and Everhart 15 have shown that the width 
of p-n junction depletion layers can be measured using the sem. 
Furthermore, Child, Ranasinghe and White 16 have shown that the 
depletion layer in mos transistors can be directly measured using the 
sem in the beam-induced conductivity (bic) mode. 1 " The MacDonald 
and Everhart technique was used on angle-lapped samples, but the 
resolution obtained was not adequate for the present study. The bic 
mode of exposure was used on laser-scribed samples and provides 
acceptable images; however, scattering of the incident electron beam 
and collection of bic electrons away from the depletion layer by 
diffusion make interpretation of the images unduly complex for junc- 
tion profile measurements. 1 ' 

The most successful method of profile measurement involves differ- 
ential etching of doped material after cleaving or angle lapping. The 
most successful etch used was nitric acid containing 1 to 2 percent 
hydrofluoric acid. Etching for 2 to 3 seconds under strong illumination 
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is adequate. A 5 percent hydrogen-peroxide solution with 1 to 2 percent 
hydrofluoric acid was also used successfully with 3-minute etch time. 
The nitric-based solution has the added advantage of acting as an 
optical stain which makes sample preparation easier. After etching, 
samples are not operated electrically, but are coated with 100 to 200 A 
of gold to improve electron emission. 

The electron microscope used is equipped with image-processing 
equipment which allows contrast ratios in the final image to be 
improved by contrast adjustments of pseudo-color image processing. 
A typical image is shown in Fig. 22. The thin rectangle at the top is 
the polysilicon gate electrode with a width of 3.43 fim ± 0.14 nm. 
Lateral diffusion of the arsenic source and drain implants under the 
gate results in the source and drain junctions etched out below the 
gate. Vertical junction depth is 0.57 u ± 0.025 /un; the ratio of lateral- 
to-vertical junction depth is 0.66 ± 0.09. These means and standard 
deviations were obtained by making several sem photographs and 
calculating the means and standard deviations from this distribution 
of measurements. 

For those lengths which are resolvable by optical means, the length 
ratios given above are within one standard deviation of the sem data. 
Optical measurements of our diffusion are not possible, however, and 




Fig. 22— sem photograph of an exposed junction. 
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in general measurements on the sem have at least five times the 
resolving power of the optical measurements, greater depth of field, 
and freedom from optical interfacing effects in lapped oxides. 

The diffusion of the arsenic implant for the samples studied was 
numerically simulated. As described in Section 5.3, the initial implant 
was modeled by a Gaussian with a peak concentration of 114.39 
(relative to the intrinsic carrier concentration) and a spread of 0.017 
/urn. The implant was then diffused for a period of 900 seconds at a 
temperature of 1373.16°K. In the samples studied, the ratio of surface 
concentration to bulk concentration was 10 4 . Applying inverse linear 
interpolation to the results of the simulation (see Figs. 10 and 14), we 
calculate that the lateral junction depth is 0.486 /im ± 0.021 /tm, the 
vertical junction depth is 0.687 urn ± 0.018 jum, and the ratio of lateral 
to vertical junction depths is 0.707 ± 0.051. The lateral junction depth 
is computed from the two-dimensional solution, while the vertical 
junction depth is computed from the more accurate one-dimensional 
solutions. Thus, the nonlinear model predicts greater penetration than 
observed experimentally. This is consistent both with the physical 
assumptions presented in Section II and with the fact that the simu- 
lation neglects oxide growth during the drive-in cycle. However, the 
shape of the diffusion front, as evidenced by the ratio of the lateral-to- 
vertical junction depths, is in excellent agreement with the predictions 
of the nonlinear model. 

VIII. PROCESS SIMULATION 

The results presented in this paper demonstrate that accurate 
simulations of realistic diffusions involving a single impurity can be 
obtained with limited computing resources. Clearly, numerical simu- 
lations should lead to deeper understanding in more complex process- 
ing situations such as multiple diffusions of interacting species and 
diffusions involving narrow masks and/or windows. Such simulations 
may, in fact, prove essential in developing technologies. 

With respect to the limited computing resources, all the simulations 
presented in this paper were performed on two minicomputers, the 
Harris Slash 7 and the Interdata 8/32. Both provide 256K words of 
directly addressable primary memory. In the case of the Harris Slash 
7, this was accomplished through a virtual memory system. In our 
implementation, which relies on direct banded-matrix techniques, 
primary memory has proven to be the crucial resource. Typical run 
times for the fine mesh solutions have been on the order of one to two 
hours. However in the case of high concentrations of donor impurities 
the run times have been on the order of 10 to 16 hours. Several months 
after the completion of this study, some of the simulations were rerun 
on a Cray-1 computer, the largest and fastest commercially available 
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scientific computer. The most difficult of the simulations required only 
7.25 minutes. Thus, in a production environment simple economics 
will favor large computers. 

Finally we note that if equilibrium solutions exist for (16), then the 
constant-surface-concentration diffusion problem for the ideal semi- 
infinite device has a family of self-similar solutions, parameterized by 
the surface concentration, « () . For a suitable selection of a ( ,'s, such 
equilibrium solutions could be calculated numerically. The resulting 
collection of solutions could then supplement the linear model pre- 
sented in Ref. 2. 

IX. CONCLUSIONS 

In this paper we investigated the combination of nonlinear and 
geometric effects on the impurity atom distribution by studying a two- 
dimensional, concentration dependent diffusion model in a region 
containing the edge of the diffusion mask. This study involved four 
major efforts. First, we developed software for solving the nonlinear 
partial differential equation describing the diffusion model, and we 
used this software to perform several numerical simulations on two 
different minicomputers, the Harris Slash 7 and the Interdata 8/32. 
Second, we applied a number of checks on the accuracy of the results 
from the simulations and concluded that the results are sufficiently 
accurate for device design using current fabrication technologies. In 
particular, we find that large minicomputers are adequate for process 
simulation of this limited scope. Third, we analyzed the simulation 
results in the light of the behavior of the nonlinear diffusion coefficient, 
D(a). From its behavior, one would anticipate that the nonlinear 
effects would become increasingly pronounced as the ratio of the 
impurity concentration to the intrinsic carrier concentration became 
larger, and that they should be more pronounced for donor impurities. 
The simulations confirmed this and demonstrated that there are three 
principal nonlinear effects: (i) a large translation of the diffusion front, 
(«) a marked steepening of the front itself, and (Hi) a very noticeable 
decrease in the ratio of the lateral-to-vertical diffusions. One inescap- 
able consequence is that, as the mask width approaches the junction 
depth, channel widths will be determined by these nonlinear, two- 
dimensional effects. Fourth, one of the authors, C. L. Wilson, obtained 
experimental measurements of actual lateral-to-vertical diffusion ra- 
tios for arsenic diffused in silicon. These measurements confirm the 
validity of the concentration dependent diffusion model. 
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APPENDIX 

In this appendix we present the details of our semi-discretization 
algorithm for solving problems such as the nonlinear diffusion equa- 
tion. We start with a partial differential equation in general divergence 
form. 

f(t, X,y, U, U X , Uy, U,, Uxt, Uyl) = D X g\{t, X,y, U, U X ,Uy,U,, U Xl , Uy,) 

+ Dyg-At, x, y, u, u x ,u y ,Ui, Uxt, u yt ) (17a) 

with initial conditions 

u(0,x,y)=h(x,y) (17b) 

and boundary conditions 

a{t,x,y)n(x,y)-(gi,g2) = b(t, x,y, u, u t ), (17c) 

where the unknown function, u(t, x, y), is defined on the domain 
ss t =s fsiop, xi =S x *s x h , and y, =s= y «s y h , and where the vector, n{x, y), 
denotes the outward pointing normal along the boundary of the 
rectangle determined by Xi, x h ,yi, and yh. 
Given meshes 

XL = Xo < X\ < • ■ • < X m -\ < X m = XR 

and 

yL = yo<y\< ■■• < y»-\ < y« " ^r» 

we will approximate u(t, x, y), the solution of (13), by a linear combi- 
nation of second-order B-splines 

u(t,x,y)= "i i U jk {t)Bj{x)Bdy). (18) 

A second-order B-spline, Bj(x), also called a chapeau function, is a 
piecewise linear function which is 1 at xj and at all other mesh points. 
Outside the interval [x L , x«], the spline is defined to be 0. Inside the 
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interval, it is defined by requiring it to be continuous. The set of 
product functions {Bj(x)B,.(y):0 ^j^m, O^k^n) is commonly 
referred to as the tensor product basis. 

The basic idea is to choose the coefficients, U j/c (t), so that the 
remainder function, 

r(t, x, y, u) = f(t, x, y, u, u x , u y , u,, u xt , Uy, ) 

- Dxg\(t t X,y,ll, Ux, Uy, Ui, Uxt, Uy,) 

- Dyg 2 (t, X, y, U, U x , Uy, U,, Uxt, Uy,), 

is orthogonal to the basis elements {Bj(x)B^{ v)}. Thus we obtain the 
following system of equations 

[/- Dxgi - D y g 2 ]Bjk dx dy = 0, (19) 

R 

where j = 0, ■ ••, m; k « 0, • • ■ , n and By* denotes the basis function, 
Bj( x)Bkiy) . This is a system of(m + l)x(n+l) nonlinear equations 
in the (m + 1) X (n + 1) unknowns, U Jk (t). We can rewrite (19) as 



[fBj, + (DxBjk)gi + iD,Bjk)gi] dx dy 



R 



D x (B Jk g x ) + Dy(B Jk gi) dx dy. 

IR 

Applying Green's theorem to the right-hand side, we obtain 
[fBj k + (DxB jk )g x + (D y B jk ) g2 \ dx dy 



R 



= n>(BjtguBjkgt) ds. (20) 



By the definition of the B j/: , the right-hand side of (19) will vanish 
whenever j < m or < k < n. However, for those splines which do 
not vanish on T, this integral will enable us to incorporate the boundary 
conditions (17c) into the Galerkin formulation. 

The next step is to replace the integrals in (19) with simple quad- 
rature formulas. For a detailed discussion of the technical details, see 
Ref. 11. This reduces (20) to an (m + 1) X (n + 1) implicit system of 
ordinary differential equations which may be written in the form 

FU, U, IT) = 0, 
with 

U(0) = U„. (21) 
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Systems of ordinary differential equations derived from such semi- 
discretization processes are usually stiff} 2 The solution of stiff equa- 
tions typically requires the use of an implicit method. Any of a number 
of standard techniques 13 can be applied. We start with the strongly A- 
stable implicit Euler rule. Then, using Schryer's step-size and order 
monitor, 18 we extrapolate trial answers obtained by this method to 
obtain a powerful variable-order, variable-step time evolution. 

This reduces the problem to the repeated solution of a large system 
of nonlinear equations. These are solved using a variant of Newton's 
method in which the Jacobian matrix is only infrequently recomputed. 
The resulting linear systems are solved with direct methods. This 
approach insures that the nonlinear iterations converge quickly at the 
expense of large storage requirements for the factored Jacobian matrix. 
The use of direct techniques also allows us to ignore certain mesh 
restrictions which might be necessary to insure the convergence of a 
fully iterative technique such as successive overrelaxation. As noted in 
Section IV, this additional robustness is particularly important during 
the initial transient analysis on the nonlinear diffusion problem. 
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